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Abstract 



We determine the second, third, and fourth virial coefficients appearing in 
the density expansion of the osmotic pressure II of a monodisperse polymer 
solution in good-solvent conditions. Using the expected large-concentration 
behavior, we extrapolate the low-density expansion outside the dilute regime, 
obtaining the osmotic pressure for any concentration in the semidilute region. 
Comparison with field-theoretical predictions and experimental data shows 
that the obtained expression is quite accurate. The error is approximately 
1-2% below the overlap concentration and rises at most to 5-10% in the limit 
of very large polymer concentrations. 

PACS: 61.25.Hq, 82.35.Lr 



I. INTRODUCTION 



For sufficiently high molecular weights, dilute and semidilute polymer solutions under 
good-solvent conditions exhibit a universal scaling behavior. 1-4 For instance, the radius of 
gyration R g , which gives the average size of the polymer, scales as N u , where N is the degree 
of polymerization and v a universal exponent, v ~ 0.5876 (Ref. 5). The osmotic pressure II 
is one of the most easily accessible quantities in polymer physics. When N is large, it obeys 
a general scaling law 2-4 (here and in the following we only consider monodisperse solutions) 

where c is the polymer number density, p the ponderal concentration, M the molar mass of 
the polymer, and T the absolute temperature. The function f(x) is universal, so that the 
determination of Z in a specific model allows one to predict II for any polymer solution. In 
the dilute limit the compressibility factor Z can be expanded in powers of c obtaining 6 

Z = l + Y J B n+ iC n 1 (1.2) 

71=1 

where the coefficients B n are knows as virial coefficients. Knowledge of B n allows one to 
compute II in the dilute regime in which cR 3 <C 1. The coefficients B n depend on the 
polymer solution. Eq. (1.2) can be rewritten as 

Z = 1 + J2 A n+i(cR 3 g ) n A n+1 = B n+1 R~ 3n , (1.3) 

71=1 

where R g is the zero-density radius of gyration of the polymer. The general scaling law (1.1) 
implies that for N — > oo the coefficients A n+ i approach universal constants A* n+1 that are 
independent of chemical details. The value of A\ has been the object of many theoretical 
studies (sometimes, the interpenetration radius \& = 2(47r)~ 3 / 2 A 2 is quoted instead of A 2 ). 
The most precise estimates have been obtained by means of Monte Carlo (MC) simulations: 
A* = 5.494 ± 0.005 (Ref. 7), A* = 5.504 ± 0.007 (Ref. 8), A* 2 = 5.490 ± 0.027 (Ref. 9). 
Eq. (1.1) is only valid for high molecular weights. Since in many cases iV is not very large, it 
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is important to consider also the leading correction to this expression. As predicted by the 
renormalization group and extensively verified numerically, for N large but finite we have 

+ A» lS ^j±„A^ l + *g£, (1.4) 

where b 2 — 1, A is a universal exponent whose best estimate is 10 A = 0.515 ± 0.0071q;qJq, 
and, of course, g(x) = 1 + J2b n+1 x n . The function g(x) as well as the constants b n are 
universal. All chemical details as well as polymer properties — for instance, the temperature — 
are encoded in a single constant a 2 that varies from one polymer solution to the other. In 
this paper we extend the previous calculations to the third and fourth virial coefficient, 
computing A3, A\, and 63. For this purpose we perform an extensive MC simulation of the 
lattice Domb- Joyce model, 11 considering walks of length varying between 100 and 8000 and 
three different penalties for the self-intersections. 

Knowledge of the first virial coefficients and of the leading scaling corrections allows us 
to obtain a precise prediction for the osmotic pressure in the dilute regime (the expression 
is apparently accurate up to B 2 c ss 1), even for relatively small values of the degree of 
polymerization. Finite-length effects are taken into account by properly tuning a single 
nonuniversal parameter. Once the virial expansion is known, we can try to resum it to 
obtain an interpolation formula that is valid in the semidilute regime. We will show that a 
simple expression that takes into account the large-density behavior of Z provides a good 
approximation to Z, even outside the dilute regime. 

The paper is organized as follows. In Sec. II we derive the virial expansion for a polymer 
solution. In Sec. Ill we define the model, while in Sec. IV we compute the universal constants 
defined above that are associated with the virial coefficients. In Sec. V we present our 
conclusions and, in particular, give an interpolation formula for Z that is also valid in the 
semidilute regime. Some technical details are presented in the Appendix. 
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II. VIRIAL EXPANSION 



We wish now to derive the virial expansion for a polymer solution. Such an expansion 
is easily derived in the grand-canonical ensemble. 12 For this purpose, we first define the 
configurational partition function Q L of L polymers in a volume V: 



Ql = J rf/iexp (^Pj2K tei ^J 



d\x 



IP 



exp 



intra 



(2.1) 
(2.2) 



where V^ ter is the sum of all terms of the Hamiltonian that correspond to interactions 
of monomers belonging to two different polymers i and j, y^ intra is the contribution due 
to interactions of monomers belonging to the same polymer i, and r ai is the position of 
monomer a belonging to polymer i. Then, we set 



(2.3) 



The virial expansion is obtained by performing an expansion in powers of z. We introduce 
the Mayer function 



f i3 = exp(-pVi^) - 1 



(2.4) 



and define the following integrals: 



'^12^13 



ri2,ri3,ri4 



h = J d 3 r 12 (/i 2 )o,r 12 
h = J d 3 r 12 d 3 r 13 (712/13/23)0 

h = J dh u dh 13 dh u (f u f 23 f 3 J u (3 + 6f 13 + f 13 f 24 )) 

r -1 2 

Ti = y rf 3 ri 2 rf 3 ri3 (/i 2 /i3) , ri2 , ri3 - [/ ^12 (/i 2 ) , ri2 

T 2 = J dh l2 dh l3 dh lA (/i 2 /i3/i4)o, ri2 , ri3 , ri4 - J d 3 r 12 (/ 12 ) 0; 

T 3 = J dh 12 d 3 r 13 dh u (/i2/23/34) , ri2 , ri3 , ri4 - J d 3 r 12 (/ 12 ) 0j 
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(2.5) 
(2.6) 
(2.7) 

(2.8) 

(2.9) 

(2.10) 



I ,ri2,ri3,ri4 



T 4 = J d 3 r 12 d 3 r 13 d 3 r u (/^a/is/s^o, 

d 3 r 12 (/i 2 ) , ri2 y d 3 r 12 d 3 r 13 (712/13/23)0, 



ri2,ri3 



(2.11) 



Here (-)o,r indicates an average over two independent polymers such that the first one starts 
at the origin and the second starts in r. Analogously (-)o,r 2 ,r a and (•)o,r 2l r 3 ,r4 refer to averages 
over three and four polymers respectively, the first one starting in the origin, the second in 
r 2 , etc. 

Then, a simple calculation gives 



v 



\n~=z + -I 2 + -(/ 3 + 37\ + 3/ 2 2 ) 

+ — (h + 4:T 2 + 12T 3 + 12T 4 + 16/f + 12/ 2 / 3 ) + 0(z 5 ). 



(2.12) 



The density is obtained by using 



d(3U 

dz 



(2.13) 



The previous equation can be inverted to obtain z in powers of c. Substituting in Eq. (2.12), 
we obtain expansion (1.2) with 



B 2 = --h, 

B, = -\h ~ \t 2 - \T, - ^T 4 + h 2 T x . 



(2.14) 
(2.15) 
(2.16) 



Note that there are additional contributions to B 3 and B 4 which are missing in simple fluids. 12 
Indeed, for a monoatomic fluid T n = 0. These terms are instead present in the polymer virial 
expansion. 

In the following we shall consider a lattice model for polymers. In this case, the previous 
expressions must be trivially modified, replacing each integral with the corresponding sum 
over all lattice points. 
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III. MODEL AND OBSERVABLES 



Since we are interested in computing the universal quantities A* n and b n , we can use any 
model that captures the basic polymer properties. For computational convenience we con- 
sider a lattice model. A polymer of length N is modelled by a random walk {r , r 1? . . . , r^} 
with |r a — r a+ i| = 1 on a cubic lattice. To each walk we associate a Boltzmann factor 

e^ H = e~ w °, a= £ 5 ra , rp , (3.1) 

0<a<f3<N 

with w > 0. The factor a counts how many self- intersections are present in the walk. This 
model is similar to the standard self-avoiding walk (SAW) model in which polymers are 
modelled by random walks in which self-intersections are forbidden. The SAW model is 
obtained for w = +oo. For finite positive w self-intersections are possible although energeti- 
cally penalized. For any positive w, this model — hereafter we will refer to it as Domb- Joyce 
(DJ) model — has the same scaling limit of the SAW model 11 and thus allows us to compute 
the universal scaling functions that are relevant for polymer solutions. The DJ model has 
been extensively studied numerically in Ref. 10. There, it was also shown that there is a 
particular value of w, e~ w * as 0.603 (i.e., w* 0.505838), for which corrections to scaling 
with exponent A vanish: the nonuniversal constant a<i is zero for w = w*. Thus, simulations 
at w = w* are particularly convenient since the scaling limit can be observed for smaller 
values of N. 

In the simulations we measure the virial coefficients using Eqs. (2.14), (2.15), and (2.16). 
In this model the Mayer function is simply 

fij = e-^r - 1 = e -*« - 1 (3.2) 
(T ij='52 S '<*i,'pr ( 3 - 3 ) 

a/3 

Here r ai is the position of monomer a of polymer %. The DJ model can be efficiently simulated 
by using the pivot algorithm. 13-16 For the SAW an efficient implementation is discussed in 
Ref. 17. The extension to the DJ model is straightforward, the changes in energy being 
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taken into account by means of a Metropolis test. Such a step should be included carefully 
in order not to loose the good scaling behavior of the CPU time per attempted move. We use 
here the implementation discussed in Ref. 18. The virial coefficients have been computed by 
using a simple generalization of the hit-or-miss algorithm discussed in Ref. 8. Some details 
are reported in the Appendix. 



IV. MONTE CARLO DETERMINATION OF THE VIRIAL COEFFICIENTS 

We perform three sets of simulations at w — 0.375, 0.505838, 0.775, using walks with 
100 < N < 8000. Results are reported in Tables I, II, and III. As expected, the results 
for w = 0.505838 are the least dependent on A, confirming that for this value of w scaling 
corrections are very small. On the other hand, for w = 0.375 and w = 0.775 scaling 
corrections are sizable. 

The numerical data are analyzed as discussed in Ref. 7. We assume that A n has an 
expansion of the form 

, ,„ T % a n (w) c n (w) 

A n (N, W )*A: + ^ + J ^L. (4.1) 

For A we use the best available estimate: A = 0.515 ± 0.007+o;oqo (Ref. 10). The term 
l/A A2 - off should take into account analytic corrections behaving as 1/N, and nonanalytic 
ones of the form N~ 2A , N~ A2 (A 2 is the next-to- leading correction-to-scaling exponent). As 
discussed in Ref. 7, one can lump all these terms into a single one with exponent A 2jC fr = 
1.0 ± 0.1. In order to estimate A 2 we also use the results that appear in Table 1 of Ref. 7 
that refer to MC simulations of interacting SAWs. 19 A combined fit of the results for A 2 
(with 1 + 5x2 = 11 free parameters) gives 

A* 2 = 5.4986 ± 0.0010 ± 0.0002 A min = 250, 

A* = 5.4997 ± 0.0017 ± 0.0002 A min = 500, 

A* 2 = 5.5013 ± 0.0030 ± 0.0003 A min = 1000 . (4.2) 
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In each fit we have considered only the data with N > iV mm . We do not show results for 
-Wmin — 100 since in this case the fit has a somewhat large x 2 - We report two error bars. 
The first one is the statistical error while the second gives the variation of the estimate as 
A and A 2je fr vary within one error bar. The results show a small upward trend which is in 
any case of the order of the statistical errors. As our final estimate we quote 

A* 2 = 5.500 ±0.003. (4.3) 

Note that in the polymer literature one often considers the interpenetration ratio ^ = 
2(47r)" 3 / 2 A 2 instead of A 2 . We have 

tf* = 2(4n)- 3/2 A* 2 = 0.24693 ± 0.00013. (4.4) 

The same analysis — but in this case we only rely on the DJ results — can be repeated for 
A 3 . Estimates of A 3 are reported in Table II. The contribution proportional to I3 appearing 
in Eq. (2.15) — the only one present in simple fluids — accounts for most of the result since 
the contribution proportional to Ti is small, Ti-R~ 6 ks 0.84 in the scaling limit. Still in a 
high-precision calculation, 7\ cannot be neglected, giving a 9% correction. A fit of the results 
to Eq. (4.1) gives 

Al = 9.788 ± 0.005 ± 0.002 7V min = 100, 

Al = 9.786 ± 0.008 ± 0.001 JV min = 250, 

Al = 9.798 ± 0.013 ± 0.001 JV min = 500, 

A* 3 = 9.813 ± 0.015 ± 0.001 JV min = 1000 , (4.5) 

where, as before, the first error is the statistical one while the second is related to the error 
on A and A c g. As final estimate we quote 

Al = 9.80 ± 0.02 . (4.6) 

In the theoretical literature, several estimates have been reported for the large-A^ value of 
^3,22 = B3/B2 (this quantity is often called g), which is universal and independent of the 
radius of gyration. Using Eqs. (4.3) and (4.6) we obtain 



^3,22 = A ll( A l) = 0-3240 ± 0.0007. (4.7) 

Finally, we consider A 4 . In this case statistical errors are quite large. This is due to significant 
cancellations among the different terms appearing in Eq. (2.16). Moreover, while in A 3 the 
term T\ was providing only a small correction, here inclusion of the terms proportional to 
Tj is crucial to obtain the correct result. They are not small: in the scaling limit we have 
T 2 R g 9 sa 28, T 3 R g 9 « 18, T A R g 9 « 1.5, TJ 2 R- 9 « -9. Fits to Eq. (4.1) give 

A\ = -9.2 ± 0.5 N mhl = 100, 

A\ = -9.3 ± 0.7 N mhl = 250, 

A\ = -9.6 ± 1.3 AUn = 500. (4.8) 

The systematic error is negligible in this case. In order to improve the result we have 
repeated the analysis taking into account that a 4 (u>) = & 4 a 2 (u>), with 6 4 independent of w. 
If we analyze together the data for A 2 and A A taking as free parameters A 2 , A* 4 , 6 4 , a 2 (w;), 
c 2 (wi), and c 4 (u>j), the nonlinear fit gives: 

A\ = -8.89 ± 0.29 AUn = 250, 

A\ = -9.00 ± 0.36 N min = 500. (4.9) 

Correspondingly, we obtain 6 4 = — 9 ± 7 and 6 4 = — 2 ± 11. Comparing all results we obtain 

A\ = -9.0 ±0.5. (4.10) 

Our result for A 2 is in good agreement with previous MC estimates: A* 2 = 5.494 ± 0.005 
(Ref. 7), A* 2 = 5.504 ±0.007 (Ref. 8), A* = 5.490 ±0.027 (Ref. 9). Direct MC calculations 19 
of A3 provided only the order of magnitude since they did not consider the contribution 
proportional to T\. Ref. 20 quotes an estimate of A 3y22) A 3y22 ps 0.40-0.45. This value is 
somewhat higher than what we obtain here, but it must be noted that very short walks were 
considered (N < 50). Thus, those results are probably affected by strong scaling corrections. 
Overall, our estimate of A 3>22 is in agreement with the field-theory estimates 21 ' 22,4 that vary 



between 0.28 and 0.43. A field-theory estimate of A\ can be obtained from the expressions 23 
presented in Sec. 17.3.2 of Ref. 4. The result, A\ ss —30, is somewhat too large, but it at 
least agrees in sign with ours. 

The fits reported above give also the coefficients a n (w) = a 2 (w)b n . In Ref. 10 it was 
claimed that a 2 (w) « for w = 0.505838. We can verify here this result. More importantly, 
we can test the renormalization-group prediction a n (w) = a 2 (w)b n , by verifying that not 
only does a 2 (w) approximately vanish, but that the same property holds for the coefficient 
az(w) [we are not precise enough to estimate reliably a^iw)}. From the fits we obtain for 
w = 0.505838: 

a 2 (w) 



A* 



0.08 ±0.02, (4.11) 



^ =0.2*0.1. (4.12) 

^3 

These results are not fully consistent with those of Ref. 10. Indeed, we find that, for w = 
0.505838, a n (w) is not zero within error bars. Still, w = 0.505838 is very close to the optimal 
value w* for which a 2 (w*) = 0. Indeed, for w = 0.505838 corrections are a factor-of-10 
smaller than those occurring for w = 0.375 and w = 0.775 and a factor-of-20 smaller that 
those occuring in SAWs with (5 — and f3 — 0.1, the values used in Ref. 7. Our best estimate 
for w* is w* = 0.48 ± 0.02. We have also recomputed the optimal p introduced in Ref. 7, 
which gives the optimal combination of SAW data corresponding to (3 — and (3 = 0.1. We 
obtain p opt = 0.52 ± 0.07. 

Finally, we compute the universal scaling-correction coefficient 63. We use the same 
method as described in Ref. 7. We define 

o ( ] _ A n (N,w 1 )-A n (N,w 2 ) 

Rn(N) = A 2 (N, Wl )-A 2 (N t w 2 y (413) 
which should scale asymptotically as 7 

P/An . , d n {w 1 ,w 2 ) , v 

R n (N) = b n + iyAeff — , (4.14) 

with A cff = 0.5±0.1. We use the three possible choices of w\ and w 2 , verifying the universality 
of the large- N behavior of R 3 (N). In Fig. 1 we report Rs(N) for the different cases. It is 
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clear that asymptotically all quantities converge to the same value, as predicted by the 
renormalization group. A fit of the data gives 

6 3 = 4.75 ±0.30, (4.15) 

where the error includes the statistical error and the systematic error due to the uncertainty 
on A cff . In principle the same analysis can be applied to At, but here errors are so large 
that no reliable estimate can be obtained. 



V. OSMOTIC PRESSURE 

The results of the previous Section allow us to determine the osmotic pressure in the 
dilute regime. Indeed, neglecting terms of order c 4 we can write 

Z=^-^l + 1.313% + 0.559$^ - 0.122^ + k$(% + 1.13$J + h^ 3 p ), (5.1) 

where we have introduced the polymer packing fraction 

AnR 3 a 4nB%N A , x 

$ p = 9 -c= 9 -—p, 5.2 

3 3 M V ; 

R 2 g is the zero-density radius of gyration, Na the Avogadro number, M the molar mass 

of the polymer, c and p the number density and the ponderal concentration respectively. 

Equivalently, we can write 

Z « 1 + X + 0.324X 2 - 0.054X 3 + k x (X 2 + b 4 , x X 3 ), (5.3) 

where X = B 2 c, thereby avoiding any reference to the radius of gyration. The constants 
&4 5 $ and b^x depend on 64 for which we have only a rough estimate. If we trust the result 
64 = — 2 ± 11 obtained in Sec. IV, we have 64^ = —0.1 ± 0.6 and b^x = 0.4 ± 1.7. The 
parameters k$ = 3a 2 /(47rA^ A ) and k x ~ 0.0392a 2 /A^ A are nonuniversal and depend on the 
degree of polymerization. 

In Fig. 2 we plot Eq. (5.1) for A;$ = (scaling limit). It is evident that Z is linear in <3> p 
only for very small values of $ p , say $ p < 0.2. For larger values of $ p inclusion of the higher- 
order terms is crucial. Note that approximations (A3) and (A4) in Fig. 2 give very close 
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predictions for Z, indicating that in the region $ p < 1 the virial expansion gives reasonably 
accurate results, with errors of the order of a few percent. Eqs. (5.1) and (5.3) depend on a 
single nonuniversal parameter that allows us to take into account the leading corrections to 
scaling due to the finite degree of polymerization. In practice, k$ and kx can be determined 
by requiring the measured value of Z at a given value of <3> p or X to agree with expression 
(5.1) or (5.3). Then, all parameters are fixed and Eqs. (5.1) and (5.3) predict Z in the whole 
dilute regime $ p < 1, X < 1. 

To extend Eq. (5.1) or Eq. (5.3) into the semidilute regime, we must modify them to take 
into account the asymptotic behavior in the scaling limit 2 

Z ~ $V(3-1) „ $1,311 „ X 1.311 (5 4) 

Moreover, a proper resummation is necessary. Since the expansion of Z 3v ~ l = Z - 763 is 
alternating in sign, we will resum this quantity by using a Pade approximant that behaves 
as <J> P for large concentrations. We write therefore (for k<$> = 0) 

/l + 1.5260$ p + 0.7954$ 2 \ 1311 

Z w £ . 5.5 

V l + 0.5245$ p J y ' 

The numerical coefficients have been obtained by requiring this expression for Z to reproduce 
Eq. (5.1) to order $ 3 . An analogous expression in terms of X can be obtained by using 
X PS 1.313$ p in the scaling limit. 

Eq. (5.5) is of course very accurate in the dilute regime since it reproduces exactly 
the virial expansion (5.1) (see Fig. 2). We must now assess the error for larger values of 
$ p . We first compare with the renormalization-group predictions of Ref. 24. They give a 
simple parametrization of their one-loop e-expansion results for the compressibility K = 
(M/RT)dTl/dp = dpil/dc, which can be measured directly in scattering experiments (Eqs. 
6 and 7 of Ref. 24 with p = 0.32 and q = 0.42). In Fig. 3 we plot ("theor") the quantity 
-^theor/-Kintcr P — 1, where i^theor is the prediction of Ref. 24 and K iTiterp the expression derived 
from Eq. (5.5). The two predictions are close and, for X — > oo, the difference is approximately 
5%. We can also compare with the field-theory results of Refs. 21, 4. For $ p — > oo, they 
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predict Z « 1.51$J- 31 (Ref. 21) and Z « 1.99$^ 309 (Ref. 4, Sect. 17.4.1), to be compared 
with Z 1.73$p 311 obtained by using Eq. (5.5). An interpolation formula 25 for Z is also 
given in Ref. 4. It is in full agreement with ours up to $ p m 3 (differences are less than 
1%). Then, the discrepancy increases, rising to 5% for <3> p ^ 10. These comparisons indicate 
that Eq. (5.5) gives a pressure that is slightly different (5-10% at most) from the field-theory 
predictions. These differences should not be taken seriously, since one-loop field-theory 
estimates have at most a 10-20% precision, as is clear from the results for A 3)22 . We can also 
compare with the numerical results of Ref. 26. Expression (5.5) describes them reasonably 
well (differences less than 10%). 

We now compare our prediction with experiments. Ref. 27 quotes Z ^ 1.50$p 32 for $ p — > 
oo, that apparently indicates that we are slightly overestimating the pressure. Note however 
that the same data give Z ^ 1 + 1.12$ p for $ p — > [compare with Eq. (5.1)], indicating 
that there are large scaling and/or polydispersity corrections. Opposite conclusions are 
reached by comparing with compressibility results for polystyrene. Ref. 24 gives an empirical 
expression for the compressibility K that fits well several sets of data for polystyrene (Eqs. 
6 and 7 of Ref. 24 with p = 0.39 and q = 0.46). In Fig. 4 we report the experimental 
results together with the prediction obtained by using Eq. (5.5) (this figure is analogous to 
Fig. 1 of Ref. 24). Our expression follows quite closely the experimental data, though the 
experimental compressibility is larger than our prediction, as can be better seen in Fig. 3 
where we report ("expt") K expt / K intCTp — 1. Again, this discrepancy should not be taken 
too seriously, since the experimental data do not satisfy the correct asymptotic behavior: 
they give Z ~ X 1 - 39 ~ $i' 39 , to be compared with the theoretical prediction (5.4). Thus, 
the discrepancy we observe could well be explained by scaling corrections and polydispersity 
effects. 

Eq. (5.1) applies of course only to situations in which the solution is in the good-solvent 
regime. Close to the 9 point, corrections are particularly strong and cannot be parametrized 
by a single coefficient In this case, one can use the strategy proposed in Ref. 7. Work in 
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this direction is in progress. 

The authors thank Tom Kennedy for providing his efficient simulation code for lattice 
self-avoiding walks. 

APPENDIX A: DETERMINATION OF THE VIRIAL COEFFICIENTS 

In order to evaluate the n-th order virial coefficient B n we need to perform a summation 
over 

Z 3(n-i)^ For thig 

purpose we use the hit-or-miss algorithm discussed in Ref. 8 for B 2 . 
The algorithm can be trivially generalized to higher-order virial coefficients. We consider a 
walk W, with monomers r ,. . ., r^, starting at the origin (r = 0), and define a^iW) and 
aj{W) as the maximum and minimum value of the j-th coordinate among the points of the 
walk. Then, given two walks W\ and W 2 we define 

5*12 = Ki2,Mi,i 2 ] x [m 2 ,i2, M 2 ,i 2 ] x [m 3) i2, M 3tl2 ] 

= [ar(^i) - <*t(W 2 ), af(W 1 ) - a^(W 2 )} x [a 2 (W 1 ) - a+(W 2 ), at(W 1 ) - a^(W 2 )} 
xK(^i) -atiWMOV^ -a^(W 2 )} , (Al) 

and, given three walks W 1: W 2) and W 3: we define 

5*12,3 = ["^1,13 + mi >32 , Mi ; i 3 + Mi ;32 ] x [m 2 ,i 3 + m 2)3 2, M 2 ,i 3 + M 2i32 ] 

x [m 3) i 3 + m 3i32 , M 3) i3 + M 3j32 ] . (A2) 

It is easy to understand the rationale behind these definitions. If walk W\ starts in the origin 
and walk W 2 is translated and starts in a lattice point that does not belong to S± 2 , then W\ 
and the translated W 2 do not intersect each other, so that the corresponding Mayer function 
j\ 2 vanishes. Analogously, if walk W\ starts in the origin and walk W 2 is translated and 
starts in a lattice point that does not belong to Si 2j3 there is no translation of W 3 such that 
the translated W 3 intersects both W\ and W 2 . This guarantees that in the calculation of the 
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virial coefficient the product /13/23 always vanishes. With these definitions the sums that 
need to be computed for J 2 , ^3, and 7 4 can be written as 

h: A2(0,ri 2 ), 

ri265"i2 

Yl Yl A2(0> r i2)/i3(0,ri3)/ 2 3(r 1 2,r 13 ), 

ri2€Si2 ri3gSi3 

E E E 

ri2£-Dl2 ri 3 eDi3 ri4£Di4 

[/i2(0, r 12 )/ 2 3(r 12 , r 13 )/34(r 13 , r 14 )/ 14 (0, r 14 )(l + / 13 (0, r 13 ) + / 24 (r 12 , r 14 )) 
+/i3(0, r 13 )/ 14 (0, r 14 )/ 23 (r 12 , r 13 )/ 24 (r 12 , r 14 )(l + / 12 (0, r 12 ) + / 34 (r 13 , r 14 )) 
+/i2(0, r 12 )/ 13 (0, r 13 )/ 24 (r 12 , r 14 )/ 34 (r 13 , r 14 )(l + / 14 (0, r 14 ) + / 23 (ri 2 , r 13 )) 
+/i2(0, ri 2 )/i 3 (0, ri 3 )/i 4 (0, r 14 )/ 23 (ri 2 , ri 3 )/ 24 (r 12 , ri 4 )/ 34 (r 13 , r 14 )], (A3) 

with 



D 12 


— 5*12,3 H 512,4, 




D 13 


— 5l 3 , 2 fl Sl 3 ,4, 




D u 


= Si4, 2 H 5*14,3. 


(A4) 



Here fij(r, s) is the Mayer function computed for walk % starting in r and walk j starting in 
s. 

Eq. (A3) shows that the computation of the virial coefficients requires the calculation 
of finite sums. They can be determined by a simple hit-or-miss procedure that provides 
an unbiased estimate. For instance, in order to compute the contribution to I 2 we extract 
randomly £ vectors G Si 2 (a — 1, • • • , £), and compute 

^X>(0,rg), 

(1=1 

V(S 12 ) = (Mi,i2 - 7711,12 + l)(M 2 ,i 2 - 7712,12 + l)(M 3) i 2 - ™3,12 + l)- (A5) 

These considerations easily generalize to higher-order coefficients. 

The other contributions 7\, T 2 , T 3 , and T 4 do not require any additional work, since 
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they factorize in products of independent sums. For instance, to determine 7\ we need to 
compute 



/i2(0,r 12 )/ 13 (0,r 13 ) 



ri2,ri3 



/i2(0,r 12 ) 

T126S12 



X] /i 3 (0,r 13 ) 

T13G513 



(A6) 



The two sums are independent and can be evaluated separately as we did for the contribution 
to I 2 . 

In the calculation of the n-th virial coefficient with the hit-or-miss method we need 
to choose a point in a 3(n — l)-dimensional lattice parallelopiped. This is done by using 
3(n — 1) random numbers. For the fourth virial coefficient, 9 random numbers are needed to 
compute each contribution. If the random number generator one is using has nonnegligible 
short-range correlations (this is the case of congruential generators, see Ref. 28), the results 
may be incorrect. Therefore, we took particular care in the choice of the random number 
generator. We considered four different random number generators: a congruential generator 
with prime modulus 

y n = mod (16807y n _!, 2 31 - 1); 
a 48-bit congruential generator 

z n = mod(311672852„_i + 1,2 48 ); 
a 32-bit shift-register generator with very long period 

tn — ^n-1029 X0Rt„_2281) 

where XOR is the exclusive-or bitwise operation; the 32-bit Parisi-Rapuano generator 29 

x n = mod (x„_ 24 + a:„_ 55 , 2 32 ) r n = x n X0Rx n _ 61 . 

In order to compute the virial coefficients we must choose one or more lattice points in a 
given three-dimensional parallelopiped. For this purpose we generate four uniform random 
numbers a, a±, a 2 , and a 3 in [0,1): 
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a = mod (t n + 2y n , 2 32 ) x 2~ 32 , 
ai = mod (r„ + 2y„ +1 , 2 32 ) x 2" 32 , 
a 2 = mod (r n+1 + 2y n+2 , 2 32 ) x 2" 32 , 

a 3 = z n x 2" 48 . (A7) 

Number a is used to determine a random permutation cr of three elements. Then, we consider 
v = (vi,v 2 , v 3 ) = (a CT (i), a CT ( 2 ), a CT (3))- A random lattice point in [mi, Mi] x [m 2 , M 2 ] x [m 3 , M 3 ] 
is (mi+ L^i(Mi— TOi + l)J,m 2 + Lv2(M 2 -m 2 + l)J,m 3 + [^3(^3-^3 + 1)])- As a check, we 
computed the virial coefficients for hard spheres using the hit-or-miss method. We obtain: 

B 2 /V = 0.499993 ± 0.000069, (A8) 
B 3 /V 2 = 0.156231 ± 0.000062, (A9) 
B 4 /V 3 = 0.035785 ± 0.000072, (A10) 

where V is the volume of the sphere. These estimates should be compared with the exact 
results 12 0.5, 0.15625, 0.035869... They are in perfect agreement. Thus, we are confident 
that our final results, that are much less precise than those reported above, are not affected 
by any bias due to the random number generator. 
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TABLES 

TABLE I. Estimates of the ratio A2 for different values of N and w. 



N = 100 
N = 250 
N = 500 
N = 1000 
N = 2000 
N = 4000 
N = 8000 



w = 0.375 
5.20511(83) 
5.31755(81) 
5.37200(93) 
5.4116(14) 
5.4365(14) 
5.4545(13) 
5.4682(14) 



w = 0.505838 
5.49049(58) 
5.50395(63) 
5.50609(68) 
5.5048(10) 
5.50610(93) 
5.5027(10) 
5.5051(10) 



w = 0.775 
5.79591(85) 
5.69318(64) 
5.63798(93) 
5.5988(14) 
5.5676(14) 
5.5466(14) 
5.5324(15) 



TABLE II. Estimates of the ratio ^3 for different values of and w. 



w = 0.375 



w = 0.505838 



w = 0.775 



N = 100 
N = 250 
N = 500 
N = 1000 
N = 2000 
N = 4000 
N = 8000 



8.4794(54) 
8.9534(59) 
9.2024(65) 
9.3724(74) 
9.4760(99) 
9.595(11) 
9.647(10) 



9.7936(41) 
9.8248(49) 
9.8295(46) 
9.8214(71) 
9.8098(65) 
9.8180(71) 
9.8062(74) 



11.3308(61) 
10.7663(48) 
10.4819(74) 
10.2789(89) 
10.1219(91) 
10.021(11) 
9.945(11) 
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TABLE III. Estimates of the ratio for different values of N and w. 





w = 0.375 


w = 0.505838 


w = 0.775 


N = 100 


-8.64(33) 


-8.64(30) 


-7.82(53) 


N = 250 


-8.86(49) 


-8.61(34) 


-7.91(41) 


N = 500 


-7.84(48) 


-9.14(37) 


-8.55(64) 


at 1 nnn 
iv = 1UUU 


— 9.24(obj 


—9.23(65) 


— 1 .92(0 I ) 


iV = 2000 


-8.30(94) 


-8.96(69) 


-8.5(1.0) 


N = 4000 


-8.8(1.0) 


-8.19(70) 


-9.6(1.0) 


N = 8000 


-8.0(1.0) 


-9.92(77) 


-10.1(1.1) 
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FIGURES 
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FIG. 1. Estimates of R 3 (N) vs l/TV - 5 for three different cases: (a) w% = 0.775, w 2 = 0.375 
= 0.775, w 2 = 0.505838; (c) w x = 0.375, w 2 = 0.505838. 



23 




FIG. 2. Compressibility factor Z vs <£ p . We report four different approximations: (A2) it 
corresponds to Eq. (5.1) truncated to order 3> p ; (A3) truncation to order (A4) truncation to 
order (interp) interpolation formula (5.5). In all cases k$ = 0. 
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FIG. 3. Comparison of the compressibility K obtained by using the interpolation formula (5.5) 
(-K"interp) and two different approximations presented in Ref. 24: one is based on one-loop field 
theory ("theor"), one is obtained by fitting different sets of data for polystyrene ("expt"). 



25 




0.1 1 — — — — 1 

0.01 0.1 1 10 100 

X 

FIG. 4. Comparison of the compressibility K obtained by using the interpolation formula (5.5) 
("interp") and by fitting different sets of data for polystyrene ("expt") (taken from Ref. 24). 
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